% Create Fresh Environment
close all
clear all
clc

L = 201;                % Length of polynomial
t=linspace(-5,5,L);
x = hann(L).*sinc(t)';

% Roots of Unity (how fine to scan unit circle)
ru = 250;
N = 2*(ru-1);

rc = 1;                     % Middle ring
rp = 0.03;                  % Ring precision



exps = (0:L-1)';            % exponential factors

i = 1;
lf_rts=[];
theta = 0;
indx = 0;
lx = 0;
while(rc > rp)
    
    r = [rc+rp rc rc-rp];
    fftx = fft(x.*(r(1).^(0:L-1)'),N);
    fftx(1:end,2) = fft(x.*(r(2).^(0:L-1)'),N);
    fftx(1:end,3) = fft(x.*(r(3).^(0:L-1)'),N);

    magx = abs([fftx(end,1) ; fftx(1:ru+1,1)]);
    magx(1:end,2) = abs([fftx(end,2) ; fftx(1:ru+1,2)]);
    magx(1:end,3) = abs([fftx(end,3) ; fftx(1:ru+1,3)]);


    k = 2:ru+1;                 % indices of middle-disc without end values
    Bk = crB(k);                % value function on the middle-ring

    % locate the local minima
    locs = (Bk < crA(k-1)) & (Bk < crA(k)) & (Bk < crA(k+1)) ...
                & (Bk < crC(k-1)) & (Bk < crC(k)) & (Bk < crC(k+1)) ...
                    & (Bk < crB(k-1)) & (Bk < crB(k+1));

    indx = find(locs);              % indices of local minima
    lx = length(indx);              % number of minima found

    
    if (lx)
        theta = -2*pi/N*(indx-1);   % convert indices into complex numbers
        lf_rts = [lf_rts; exp(1i*theta)];
    end
    rc = rc - rp;
    r = [rc+rp rc rc-rp];

end
actual_rts = roots(flipud(x)); % compute roots with companion matrix
plot(lf_rts,'*'), hold on % plot LF approximations,
plot(conj(lf_rts),'*') % companion matrix roots
plot(actual_rts,'or') % and the three concentric discs
uc = exp(1i*linspace(0,2*pi,200)');
plot(uc*r,'k'), hold off
axis equal, axis off